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Abstract 


The singular value decomposition has been extensively used 
for the analysis of the kinematic and dynamic characteristics 
of robotic manipulators. Due to a reputation for being nu¬ 
merically expensive to compute, however, it has not been 
used for real-time applications. This work illustrates a for¬ 
mulation for the singular value decomposition that takes 
advantage of the nature of robotics matrix calculations to ob¬ 
tain a computationally feasible algorithm. Several applica¬ 
tions, including the control of redundant manipulators and 
the optimization of dexterity, are discussed. A detailed illus¬ 
tration of the use of the singular value decomposition to deal 
with the general problem of singularities is also presented. 


1. Introduction 


In recent years the singular value decomposition 
(SVD) has become a popular tool for analyzing the 
kinematic and dynamic properties of robotic manipu¬ 
lators (Yoshikawa 1985a; Yoshikawa 1985b). It plays 
a particularly prominent role with regard to redundant 
manipulators, both in terms of analyzing the signifi¬ 
cance of the extra degrees of freedom (Klein and 
Huang 1983) and in specifying a side criterion that 
can be optimized using these redundant degrees of 
freedom. In many cases, these side criteria are some 
quantitative measure of the qualitative concept of 
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dexterity. Most of the dexterity measures proposed are 
some function of the singular values of the Jacobian 
matrix. The most common of these is perhaps the 
manipulability measure proposed by Yoshikawa 
(1984) that is defined as the square root of the deter¬ 
minant of the matrix JJ T , which is simply the product 
of the singular values of /. Other proposed measures 
include the trace of the above matrix (Baillieul 1987), 
the minimum singular value of the Jacobian (Klein 
and Blaho 1987), the compatibility index (Chiu 1987), 
and isotropy (all equal singular values; Salisbury and 
Craig 1982). 

While all of the above measures have a physical 
significance and justification for their use, the key 
point here is that they are all closely linked to the 
SVD. Yet in spite of this fact, the full decomposition 
is usually limited to the analysis of manipulator con¬ 
figurations and is not considered for implementation 
in on-line control. This is exemplified by the popular¬ 
ity of the manipulability measure, since its major jus¬ 
tifications are that it is numerically simple to compute 
and that its zeros coincide with the singularities of the 
Jacobian. The implication is that one would really like 
information about singularities; however, that would 
require calculation of the SVD, which has a reputation 
for being numerically expensive to compute. Unfortu¬ 
nately, the determinant gives no information about 
the absolute proximity to singularities, since the mini¬ 
mum singular value is the only reliable measure of 
this quantity. In addition, the calculation of the matrix 
product JJ T squares the condition number, which 
reduces the accuracy of the result. 

This work is concerned with demonstrating that, 
with the right formulation, the SVD is computation¬ 
ally feasible for use in real-time control. Traditionally, 
the computation of the SVD of an arbitrary matrix is 
an iterative procedure, so that the exact number of 
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computations cannot be known a priori. However, the 
control of robotic systems is not based on the solution 
of arbitrary matrix equations but quite frequently 
involves the solution of equations based on the Jaco¬ 
bian matrix. The current Jacobian for a system can be 
regarded as a perturbation of a previously known ma¬ 
trix for which perturbation bounds on the singular 
values and singular vectors can be established. It will 
be shown how this knowledge of the previous state can 
be exploited during the current computation of the 
SVD in order to reduce the overall computational 
burden. The above point will be emphasized in this 
work, which illustrates a computational scheme capa¬ 
ble of calculating the SVD of the Jacobian for use in 
real-time control of manipulators. The implications of 
such an algorithm for several applications, including 
the utilization of redundancy and the optimization of 
manipulator dexterity, are then discussed. A detailed 
examination of the advantages of using the SVD for 
dealing with manipulator singularities is illustrated 
through computer simulations of a PUMA robot, 
using both rate and acceleration control. 


2. Overview of SVD Algorithms 


2.1. The Golub-Reinsch Algorithm 


The most popular technique for computing the SVD 
was first proposed in Golub and Kahan (1965) with 
source code published in Businger and Golub (1969). 
It is now commonly referred to as the Golub-Reinsch 
algorithm (Golub and Reinsch 1970) and is available 
in many linear algebra software packages such as 
EISPACK, LINPACK, and IMSL. This algorithm is 
composed of two distinct steps, namely, transforming 
the given matrix into bidiagonal form using a series of 
Householder transformations and then applying an 
iterative procedure designed to use orthogonal trans¬ 
formations to produce bidiagonal matrices that are 
successively more diagonal. The procedure used is a 
variant of the QR algorithm (Watkins 1982), with 
origin shifts designed to improve the convergence 
properties (Stewart 1970). 


The Golub-Reinsch algorithm is generally regarded 
as the most efficient and numerically stable technique 
for computing the SVD of an arbitrary matrix. How¬ 
ever, there are two aspects of the algorithm that make 
it less attractive for the problem at hand. The first 
relates to the fact that the first step in the algorithm 
requires a fixed computation to bidiagonalize the given 
matrix so that a slightly perturbed matrix must still 
undergo this operation. The second aspect relates to 
the relatively serial nature of the technique, thus mak¬ 
ing it difficult to utilize parallel computing structures. 
For these reasons, the Golub-Reinsch algorithm will 
not be considered as the basis for implementing a 
real-time SVD algorithm. The following algorithm, 
based exclusively on Givens rotations, is more suited 
to take advantage of incremental perturbations and 
parallel architectures. 


2.2. Algorithms Based on Givens Rotations 


Givens rotations are orthogonal transformations of the 
form 



' J 


where all other elements not shown are zero. This 
transformation can be geometrically interpreted as a 
plane rotation of d in the i-j plane. These transforma¬ 
tions are also known as Jacobi rotations since they 
were first described by Jacobi (1846). In contrast to 
Householder reflections, Givens rotations only affect 
two rows or columns of the matrix with which they 
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are multiplied, a property that is useful when designing 
parallel computing structures based on this transfor¬ 
mation. One popular application of Givens rotations is 
for the computation of the QR decomposition. While 
this requires 2mn 2 — 2/a 3 /3 floating-point operations 
as compared to mn 2 — n 2 /3 floating point operations 
for the Householder version of the QR decomposition 
(Golub and Van Loan 1983), the Givens rotations can 
be done in parallel and can therefore be computed in 
O(m) units of time with an appropriately connected 
array of 0(n 2 ) processors (Luk 1986a). 

The most important property of Givens rotations 
for the problem at hand is their ability to orthogonalize 
the two rows or columns on which they operate. This 
forms the basis of an SVD algorithm that relies solely 
on Givens rotations (Hestenes 1958; Nash 1975). In 
particular, consider an orthogonal matrix V, composed 
of successive Givens rotations, such that 

AV=B (2) 

where the columns of B are orthogonal. If the columns 
of B are orthogonal, then it can be written as the prod¬ 
uct of an orthogonal matrix U and a diagonal matrix D 

B=UD 


each of which is designed to orthogonalize two col¬ 
umns. Considering the current z'th and yth columns of 
A, multiplication by a Givens rotation results in the 
new columns, a,' and nj given by 

= a, cos (0) + a, sin (0) (7) 

a/ = aj cos (0) — a t sin (0). (8) 

The constraint that these columns be orthogonal re¬ 
sults in 

a i T aj = 0 = a/ajcos 2 (0) - sin 2 (0)] 

+ (a^aj — a?a t ) sin (0) cos (0). ' ’ 

The terms in the Givens rotation matrix to achieve 
orthogonality can be computed by using the formulas 
given in Nash (1979), which are based on the quantities 


£ 

08 

II 

(10) 

q - aft, - afo 

(11) 

v = \l4p 2 + q 2 

(12) 


(3) so that for q s 0 


by letting the columns of U be equal to normalized 
versions of the columns of B, 


cos (0) — 


and sin (0) = 


V cos (0) 


and defining the diagonal elements of D to be equal to 
the norm of the columns of B 

<4 = INI- (5) 

By substituting (3) into (2) and solving for A, one obtains 



A = UDV T (6) 

which is the SVD of A. 

The critical step in the above procedure for calculat¬ 
ing the SVD is determining the orthogonal mtrix V 
that will orthogonalize the columns of A. This matrix 
is usually formed as a product of Givens rotations. 


The two sets of formulas are given so that ill-condi¬ 
tioned equations resulting from the subtraction of 
nearly equal numbers can always be avoided. 

The preceding discussion shows how to determine a 
single Givens rotation that will orthogonalize two 
columns of a given matrix. It still needs to be shown 
how the matrix V can be computed from these ele- 
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mentary rotations. If the Givens rotation to orthogo- 
nalize columns / and j is denoted by V tJ , then the 
product of a set of n(n - l)/2 rotations denoted by 

^-n(n^) 06 ) 

is referred to as a sweep (Golub and Van Loan 1983). 
Unfortunately, a single sweep will not, in general, 
orthogonalize all of the columns of a matrix, since 
subsequent rotations can destroy the orthogonality 
produced by previous ones. However, the procedure 
can be shown to converge (Nash 1975) so that Fean 
be obtained from 



(17) 


where the number of sweeps / is not known a priori. 
Convergence of the algorithm is based on completing 
an entire sweep with all of the columns being orthogo¬ 
nal. Orthogonality is measured by the parameter a. 
defined as 


_ (a? - *,) 

(a?'a l )(a J T a J ) 


dropping below a preset threshold. If for two columns 
a * s below the threshold, then the rotation is not per¬ 
formed. 

The above algorithm, by virtue of being composed 
exclusively of Givens rotations, can be highly parallel¬ 
ized, a task that has already been performed for imple¬ 
mentation on the ILLIAC IV (Luk 1980). Architec¬ 
tures specifically designed for this algorithm have also 
been proposed (Luk 1986b, Schimmel and Luk 1986) 
and can operate at about 2 fn 2 units of time per sweep 
for annXn matrix where/is the time required for a 
floating-point multiply and add. It has been shown 
that the number of sweeps required in the above algo¬ 
rithm is approximately log 2 n. In the following sec¬ 
tions it will be shown how perturbation bounds on the 
singular values and vectors of the Jacobian can be 
incorporated into this algorithm in order to reduce the 
number of sweeps required as well as the computa¬ 
tional complexity of each sweep. 


2.3. Pertubation Bounds on the Singular Value 
Decomposition 


The algorithm for computing the SVD using Givens 
rotations outlined above uses successive sweeps in 
order to make the columns of a matrix more orthogo¬ 
nal. The more orthogonal the columns are to begin 
with, the fewer the number of sweeps required for 
convergence. If one considers the current manipulator 
Jacobian to be a pertubation of the previous Jacobian 

J(t + At) = J(t) + AJ(t), (19) 

the SVD of which is known and given by 

J(t) = U(t)D(t)V r (t), (20) 

then the matrix J(t + A/)F(/) will have nearly ortho¬ 
gonal columns, provided the perturbation A J(t) is 
small relative to /(/). The foundation of the above lies 
in the fundamentally well-behaved nature of the SVD 
of a matrix. The perturbation bounds on singular 
values are very well known and easy to show (For¬ 
sythe, Malcolm, and Moler 1977): 

+ At))- Oj(J(t ))| < || AJ(t)|| (21) 

The perturbation bounds on the rotation of subspaces 
defined by singular vectors are not as widely known 
but are also well behaved (Davis and Kahan 1970- 
Wedin 1972). 


3. Implementation of a Real-Time 
SVD Algorithm 


The implementation of the SVD algorithm using 
Givens rotations and the subsequent refinements were 
all done in PASCAL on a VAX 785. Substantial test¬ 
ing for a variety of trajectories using a simulation of 
the PUMA robot have been conducted. Three repre¬ 
sentative examples are given in Figures 1, 2, and 3, 
illustrating the starting configuration of the PUMA 
robot, the desired end-effector trajectory, and the sin- 
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Fig. 1. Initial configuration 
and desired end-effector 
positions for trajectory A, 
along with the singular 
values of the Jacobian for 
each point along this trajec¬ 
tory. 




gular values of the Jacobian computed for each point 
along the trajectory. As can be seen from the singular 
value plots, these examples all pass near singularities 
at the midpoint of their trajectories. These singularities 
are the well-known wrist, shoulder, and elbow singu¬ 
larities that occur in trajectories A, B, and C, respec¬ 
tively. Trajectory C is unique in that it approaches a 
triple singularity at its midpoint. These examples have 
been chosen to illustrate the advantages of having the 
SVD available during the control of a manipulator, 
since conventional algorithms provide unsatisfactory 
performance near singular configurations. 


Fig. 2. Initial configuration 
and desired end-effector 
positions for trajectory B 
along with the singular 
values of the Jacobian for 
each point along this trajec¬ 
tory. 




In order to illustrate the computational require¬ 
ments of the basic algorithm, which does not use pre¬ 
vious information, and to provide a basis for compari¬ 
son with the modified algorithm, data on the number 
of sweeps and plane rotations required to reach con¬ 
vergence is presented in Table 1. Figure 4 shows a plot 
of these quantities for trajectory B as a representative 
example. The number of sweeps is the actual number 
of sweeps in which rotations are performed and does 
not include the final sweep in which all the columns 
are checked and determined to be orthogonal. The 
maximum number of rotations per sweep is 15 since 
the Jacobian in this case is a 6 X 6 matrix. Note that 
the computational expense is fairly uniform over var¬ 
ious configurations of the manipulator, with the slight 
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Fig. 3. Initial configuration 
and desired end-effector 
positions for trajectory C 
along with the singular 
values of the Jacobian for 
each point along this trajec¬ 
tory. 




variation in the total number of rotations primarily 
introduced by the final sweep before convergence. 

3.1. Incorporating the Previous SVD 

The above experimental data backed by the analytical 
results pertaining to the perturbation bounds on the 


Table 1. Summary of the Computational Requirements 
for Computing the SVD Using Givens Rotations and 
No Previous Information 


Trajectory 

Rotations 


Sweeps 


Min. 

Max. 

Avg. 

Min. 

Max. 

Avg. 

A * 

39 

52 

45.7 

3 

4 

3.63 

B 

35 

56 

44.9 

3 

4 

3.24 

C 

32 

45 

40.2 

3 

3 

3.00 

Total 

32 

56 

43.6 

3 

4 

3.29 


rotation of singular vectors seems to indicate that the 
majority of the computation involved in calculating 
the SVD is redundant over a given trajectory. Each 
time the SVD is computed, the columns are orthogo- 
nahzed by building the matrix V from scratch If, on 
the other hand, the value of V(t + At) i s initialized to 
HO, then a substantial portion of the work required 
to compute the current SVD can be eliminated. The 
results of using this past information in the SVD cal¬ 
culation are presented in Table 2, which illustrates the 
dramatic reduction in computational expense. The 
number of rotations is down by about a factor of three 
to approximately 15, with convergence obtained in 
virtually one sweep. The nearly uniform requirement 
of one sweep for convergence suggests removing the 
iterative nature of the algorithm by fixing the number 
of sweeps at one. This provides the additional compu¬ 
tational advantage of removing convergence tests. 

There will still exist cases where a single sweep will 
not result in convergence, thus introducing error into 
the calculated SVD. A graph of this error for SVD 
calculations along trajectory B is presented in Figure 
5. Only data for trajectory B is plotted, since the SVf 
calculations along both trajectories A and C only re 
quire a maximum of one sweep for convergence. Sc. 

error measures h ave been computed in order to 
differentiate the type of error and its source The sin¬ 
gular value error, denoted here by is a measure of 
the error between the calculated singular values, a 
and the actual singular values, a a (computed usi’ng^the 
Golub-Reinsch algorithm in the IMSL package), de- 
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Fig. 4. The number of sweeps 
and rotations required to 
compute the SVD of the Ja¬ 
cobian along trajectory B. 


SVD Computations 



fined by the equation 


= 2 K - a cf- 

I-1 


( 22 ) 


The error in the input and output singular vectors, 
denoted by V err and U„ r , respectively, is computed 
using the equations 


F OT = max(||/— VV T \\, ||/— V T V ||) (23) 

and 

U err = max (||/ - UU T ||, ||/ - U T U\\) (24) 


where the spectral norm is used (Lawson and Hanson 
1974). Finally, the error in the Jacobian, denoted by 
J err , is computed using 


m ’ 


(25) 


which is the normalized error in using the computed 
SVD as a representation for the Jacobian. 


Table 2. Computational Requirements for Computing 
the SVD Using the Previous Estimate of the Matrix V 


Trajectory 


Rotations 


Sweeps 


Min. 

Max. 

Avg. 

Min. 

Max. 

Avg. 

A 

14 

15 

4.975 

1 

1 

1.00 

B 

14 

17 

14.970 

1 

2 

1.07 

C 

13 

15 

14.831 

1 

1 

1.00 


The first point to note about the error terms plotted 
in Figure 5 is that they are all very small in magni¬ 
tude, with a maximum on the order of 0.01%. The 
second important characteristic is the fundamental 
difference between the monotonically increasing error 
of the input singular vectors, singular values, and the 
Jacobian, as compared to the error plot of the output 
singular vectors. The monotonically increasing error in 
the input singular vectors is a result of compounding 
the roundoff error of previous computations by initial¬ 
izing V to its value from the previous computation 
interval. This error in V is in turn responsible for the 
error in the singular values and the Jacobian. This 
error, however, is not carried over into the output sin¬ 
gular vectors due to the fundamental difference in the 
way that they are computed. While V is computed as a 
product of successive plane rotations, U is computed 
by normalization of the orthogonal columns of B (see 
e Q- (4)). Therefore, since the compounded error in V 
does not affect the orthogonality of the columns of B, 
the error in U is not monotonically increasing, but 
results from the additional sweep that would be re¬ 
quired for convergence. Analysis of the sweep data 
shows that the two peaks in the output singular vector 
error correspond to those Jacobians that required two 
sweeps in order to orthogonalize their columns. It is 
important to note, however, that this error is subse¬ 
quently reduced to its previous small value. 

The monotonically increasing error due to using the 
previous V matrix is still of some concern. While the 
magnitude of this error is small for this trajectory, it 
will grow without bound. Note that for repetitive tasks 
in conservative systems this is not a problem; when 
the manipulator returns to its initial starting configu¬ 
ration, the true SVD is assumed to be known, and V 
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Fig. 5. Error in computation 
of the SVD introduced from 
incorporating the previous V 
matrix and using a single 
sweep. 


Fig. 6. The singular values 
of the Jacobian for each 
point along trajectory B 
along with the maximum 
magnitude of the plane rota¬ 
tion required during the 
computation of the SVD. 


SVD Error 


Singular Values 



Trajectory B 


0 . 01 % 




0.0 

| 0 . 7 | 


can be reset. For arbitrary open trajectories, however, 
this will not be the case. A simple solution would be to 
periodically recompute the SVD with V reset to J; 
however, this will take between three and four sweeps 
and result in a nonuniform computation time inter¬ 
val. Fortunately, there is an alternative technique that 
takes advantage of the contrast between the error in 
computing V and U. Since U is not corrupted by com¬ 
pounded error, it can be used to reset V by carrying 
out the plane rotations on the rows instead of the col¬ 
umns of J. Thus if U T J = # where the rows of# are 
orthogonal, then # — DV T where the rows of V T are 
normalized rows of #, so that once again J= UDV T . 
Since V has now been computed by the normalization 


of orthogonal rows, it is not affected by the com¬ 
pounded error in the previous V. By alternating this 
procedure each interval, applying plane rotations first 
from the right on the columns of J and the next time 
from the left on the rows, one can effectively eliminate 
the buildup of error along the trajectory. A plot of the 
SVD error terms for trajectory B using this technique 
is given in Figure 5. As expected, the monotonically 
increasing error in the input singular vectors has been 
eliminated and is therefore also no longer reflected in 
the computed singular values or in J m . The error 
terms are now all of the same form, with peaks at those 
configurations where two sweeps would be required to 
orthogonalize the rows or columns of the Jacobian. 
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Fig. 7. A graph of the plane 
rotation angle required to 
orthogonalize two vectors 
versus the relative length of 
the vectors. This function is 
plotted for various values of 
the angle between the two 
vectors. 


3.2. Small Angle Approximations 


A further reduction in the computational complexity 
of the equations for computing the SVD can be gained 
by examining the nature of the plane rotations re¬ 
quired. A plot of the maximum rotation angle required 
during the computation of the SVD of the Jacobian 
for each point along trajectory B is presented in Figure 
6. Note that for the majority of the trajectory the max¬ 
imum rotation angle is very small; however, there 
exist three peaks, two of which are very sharp. From 
comparing the position of these peaks to the spacing of 
the singular values also presented in Figure 6, it is 
clear that very large rotations are required when there 
is a crossing of adjacent singular values. By comparing 
the position of the second two peaks in Figure 6 with 
those in the error plots of the singular vectors in Figure 
5, one can see that they coincide. Thus the large angles 
in the plane rotations result in an extra sweep being 
required for the algorithm to converge, thus introduc¬ 
ing the error. It may at first seem curious that the 
largest of the peaks, the first, does not produce any 
peak in the error plot. The reason for this apparent 
anomaly is that a large rotation angle is a necessary but 
not a sufficient condition to require an additional 
sweep. From examining the spacing of the singular 
values at the position of the first peak, one can see that 
the two equal singular values, ct 4 and a 5 , are separated 
from the remaining singular values. This implies that 
the large rotation required is restricted to the plane 
defined by their associated singular vectors and can 
therefore be successfully completed within a single 
sweep. In contrast, the final two peaks that result from 
the crossing of singular values a 3 with <r 4 and o i with 
cr 6 , respectively, occur at a point where these four 
singular values are closely spaced. This results in a 
greater interaction between plane rotations, therefore 
creating the necessity of an additional sweep for con¬ 
vergence and introducing error into the single sweep 
approximation. 

Since the maximum plane rotations required during 
the algorithm are very small in magnitude outside of a 
few isolated peaks, a small angle approximation is 
useful in reducing the computational effort required in 
computing the SVD. The plane rotation required to 
orthogonalize two vectors a, and a y is obtained by 


Rotation Angie to Orthogonalize Vectors 



satisfying eq. (9). By applying the double angle for¬ 
mulas, the above equation can be solved for 9 resulting 
in 


9 = 'h tan 1 


2 a ?~ a i 

afa-aft,' 


( 26 ) 


Dividing both the numerator and the denominator of 
the right hand side of eq. (26) by ||aj|| 2 results in the 
still exact equation 


9 = 


‘A tan 



(27) 


where </> is the angle between a, and aj. This form of 
the equation makes explicit the dependence of the 
plane rotation angle on both the non-orthogonality 
of the two vectors involved as well as their relative 
length. A plot of the required plane rotation to ortho¬ 
gonalize two vectors versus their relative lengths for 
various values of cos (<£) is given in Figure 7. 

One interesting point about the graph in Figure 7 is 
that when the two vectors are of equal length they will 
be orthogonalized by a rotation of 45 degrees regard¬ 
less of the angle between them. The most important 
point in terms of using a small angle approximation, 
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Fig. 8. Error in computation 
of the SVD introduced from 
using a small angle approxi¬ 
mation. 


however, is that for nearly orthogonal vectors the rota¬ 
tion angle will be small, except when there are nearly 
equal singular values. Using the previous estimate of 
the SVD to pre-orthogonalize the current Jacobian 
guarantees that this will be true. Therefore, if the two 
vectors a, and aj are not equal in length, then 0 will be 
small. For small 0 the approximations 

cos (8) ~ 1 and sin ( 6) ~ 8 (28) 

are valid. The advantages of using this approximation 
are two-fold. First, the solution of eq. (26) is vastly 
simplified to 


8 


a J a i ~ a * a i 


(29) 


so that the expensive computations of eqs. (12)-(14) 
are no longer required. Second, the calculations to 
compute the results of this plane rotation previously 
required two floating point multiplies and one addition 
per element. By using this small angle approximation, 
eqs. (7) and (8) now become 

ai = a, + a^ and af = aj — a,0 (30) 

so that the number of floating point multiplies re¬ 
quired has been cut in half. One must still consider, 
however, the case where the two vectors are of nearly 
equal length (i.e., nearly equal singular values). As 
discussed above, under these circumstances not only is 
the small angle approximation not valid, but the angle 
of rotation is at its maximum value of 45°. Fortu¬ 
nately, the cosine and the sine of 45° are equal, so that 
the reduction in floating point multiplies can still be 
achieved. Therefore, if ||a,j| ~ ||aj|| then 

cos (6) = sin (8) — V2/2 (31) 


and eqs. (7) and (8) become 


ai'=Y (ai + »j) 


and a/ = y (aj-a,), 


(32) 


which still only require one floating point multiply 
and addition per element. 

Simulation results showing the error in the com- 


SVD Error 



1 °/o 


0 

1 % 


0 

1 % 


0 

1 % 


0 


Trajectory B 


puted SVD using the small angle approximations dis¬ 
cussed above are presented for trajectory B in Figure 
8. Trajectory B is presented as an example, since it 
represents the worst case of the three trajectories. The 
form of the error terms now more closely reflects the 
plot of Figure 6, since approximations for large rota¬ 
tions will always introduce error regardless of whether 
they are restricted to a single plane. The peak errors 
are all well within 1%, with nominal values at around 
0.001%. The peak values simply reflect the inherently 
ill-conditioned nature of trying to define singular vec¬ 
tors for nearly equal singular values. The small angle 
approximation primarily introduces error to the norm 
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of the rows and columns of the rotation matrix but 
not to their orthogonality. Thus a renormalization of 
the resultant rotation matrix significantly reduces the 
errors introduced by the small angle approximation. 
While this degree of accuracy should be sufficient for 
most applications, it is important to note that manipu¬ 
lator configurations that produce peaks in the error 
terms are easily identified by examining the spacing of 
the singular values. Therefore, the small angle approx¬ 
imations can be abandoned when there are closely 
spaced singular values, thus obtaining an even higher 
degree of accuracy while retaining the computational 
advantages of the approximation throughout the rest 
of the trajectory. 


3.3. Numerical Evaluation 

A test was performed to evaluate the actual CPU time 
needed to compute the SVD of the Jacobian matrix 
for a six-degree-of-freedom manipulator. The algo¬ 
rithm using small angle approximations and the pre¬ 
vious estimate of the singular vectors, coded in 
PASCAL and executed on a VAX785 computer, re¬ 
quired 5.13 ms. In comparison, the execution time for 
the unmodified algorithm (without small angle ap¬ 
proximations or the use of previous estimates) required 
an average of 27.2 ms for a typical Jacobian. This 
figure is comparable to the time required by the 
Golub-Reinsch algorithm in the IMSL package 
(32.3 ms). These figures, however, do not reflect the 
advantage of the algorithm in terms of its parallelism, 
since they are coded and executed on a serial machine. 
By using a simple mesh connection of processing ele¬ 
ments capable of executing an addition or multiplica¬ 
tion, the plane rotations can be computed in parallel. 
Such architectures suggest that a computation time of 
well within 1 ms are easily achievable (Luk 1986b; 
Schimmel and Luk 1986). 


4. Applications 

The ability to calculate the SVD of the Jacobian in 
real time has a number of possible different applica¬ 


tions. As mentioned above, several proposed dexterity 
measures are closely linked to the singular values and 
vectors of the Jacobian. Thus the kinematic and static 
force capabilities of the current manipulator configu¬ 
ration can be compared to the requirements of an 
assigned task. Inconsistencies between the physical ca¬ 
pabilities of the manipulator and the assigned task can 
be addressed by the manipulator itself by determining 
a more suitable configuration. This results in more 
autonomous behavior, which can be utilized for both 
automated motion planning and work-cell design as 
well as for the operation of robotic manipulators in 
unstructured environments. 

Real-time computation of the SVD also enhances 
the utilization of redundancy in robotic systems. In 
terms of the resolved motion rate control formulation 
(Whitney 1969), 

J6 = x (33) 

one of the most common techniques for using the 
redundant degrees of freedom within the system is to 
use the projection operator formulation proposed in 
Liegeois (1977) 

9 = J + x + (I-J+J) z ( 34 ) 

where J + is the pseudoinverse of J and z is an arbitrary 
vector in 9 space. This formulation has been used to 
optimize a number of secondary criteria under the 
constraint of a specified end-effector trajectory includ¬ 
ing joint availability, torque minimization (Holler- 
bach and Suh 1987), and obstacle avoidance (Macie- 
jewski and Klein 1985; Nakamura, Hanafusa, and 
Yoshikawa 1987). With the complete SVD available, 
the projection operation becomes trivial, since the 
singular vectors v, for r > i 2 * n specify an orthonormal 
basis for the null space. Thus the relative advantages 
of using the homogeneous solution for alternate sec¬ 
ondary criteria can be easily evaluated. 

The remainder of this work will consider the appli¬ 
cation of the real-time SVD algorithm to the funda¬ 
mental problem of singular configurations. With the 
exception of Cartesian positioning robots, all articu¬ 
lated manipulators can be shown to possess singular 
configurations that limit the effective number of inde¬ 
pendent degrees of freedom (Baker and Wampler 
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1987). A considerable amount of effort (Asada and 
Cro Granito 1985; Aboaf and Paul 1987; Dubey and 
Luh 1987; Mayorga and Wong 1987; Sampei and 
Furuta 1987) has been devoted to either avoiding or 
dealing with operations at singularities due to the high 
joint velocities and spurious motions that can result. 


4.1. Damped Least-Squares Solutions 


The effects of singularities are frequently presented 
with respect to the resolved motion rate control for¬ 
mulation given in eq. (33). Singularities are identified 
by a mathematical change of rank in /, which physi¬ 
cally represents the inability of the manipulator to 
achieve an arbitrary end-effector velocity. For these 
cases, inverses are not defined, and even pseudoinverse 
solutions such as eq. (34) are unsatisfactory since there 
is an undesirable discontinuity at the singularity that 
can result in oscillations and unacceptably high joint 
velocities. These difficulties are not unique to the re¬ 
solved motion rate formulation but are an inherent 
part of the transformation between Cartesian and joint 
space. 

A general approach to resolving the discontinuity 
at singular configurations and maintaining a well- 
conditioned formulation that results in physically 
meaningful joint velocities is to use the damped least- 
squares formulation independently proposed in Naka¬ 
mura and Hanafusa (1986) and Wampler (1986). The 
damped least-squares solution of eq. (33) is the solu¬ 
tion that minimizes the sum ||x — J6\\ + X\\6\\ so that 
the end-effector tracking error is weighted against the 
norm of the joint velocity by using X, also known as 
the damping factor. This solution is typically obtained 
by solving an equation of the form 

(J T J+X 2 I)0 = J T x. (35) 

This solution is guaranteed to be the minimal residual 
solution over all solutions of equal or smaller norms. 
Unfortunately, the norm of the solution cannot be 
determined a priori for a given damping factor. 

The specification of the problem one would like to 
solve is the minimization of the residual ||x — J6 1|, 


which defines the end-effector tracking accuracy under 
the constraint ||0|| ^ 6 max where 9 max is the physical 
limit on the manipulator’s joint velocity. The desired 
solution can therefore be obtained by using the 
damped least-squares solution for an appropriate value 
of the damping factor. Intuitively, if the value of ||0|| 
for which the residual is equal to zero is less than dm^, 
then X = 0; otherwise X would take on the value that 
results in ||0|| = d max . In physical terms, if the joint 
velocity that exactly tracks the desired end-effector tra¬ 
jectory is physically achievable, then it should be used; 
otherwise the optimal solution requires that the joint 
velocity norm be at its limit. 

The damped least-squares solution of eq. (33), which 
will be denoted by 0 (A) in order to denote its explicit 
dependence on the damping factor, is given by 


where 




J = 2 <7 < U l V i r 

/-I 


(36) 


(37) 


is the SVD of the Jacobian. The solution norm is, 
therefore, given by 


where 



X, = u^x. 


(38) 


(39) 


Equation (38) is the nonlinear equation in A that must 
be solved in order to find the optimal solution when 
Ilia’ll = 9 max . An efficient technique for doing so is to 
use Newton’s method, which requires the derivative of 
eq. (38) with respect to the damping factor, as dis¬ 
cussed in Lawson and Hanson (1974). 

The above technique was implemented in a simu¬ 
lation of the PUMA robot for the three trajectories 
presented in Figures 1, 2, and 3. As mentioned pre¬ 
viously, these trajectories are chosen to go through 
singular configuration in order to illustrate the proper¬ 
ties of the damped least-squares solution. The results 
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Fig. 9. The norm of the 
end-effector velocity tracking 
error and the joint angle 
velocity for trajectory A and 
C using the damped least- 
squares solution. 


Damped Least Squares 



Trajectory C 



Trajectory A 


of the simulation, both the end-effector velocity track¬ 
ing error and joint angle velocity norm, are plotted for 
the three trajectories in Figures 9 and 10. The maxi¬ 
mum joint velocity norm 9^ was set at 0.009 radians 
per computation interval for all three trajectories. The 
end-effector tracking error is zero at all points along 
the trajectory, except where the Jacobian becomes 
nearly singular and the desired end-effector velocity 
has a component in the direction of the lost degrees of 
freedom. At these points the characteristic jump in the 
joint velocity is observed but is effectively clamped at 


the maximum allowable value. This prevents the spur¬ 
ious motions that are typical of manipulators passing 
through singular configurations using other methods of 
inverse kinematics. Thus the solutions are physically 
meaningful even near singular configurations, and, 
in fact, they result in the minimum amount of end- 
effector tracking error. This error can actually be zero 
if the commanded end-effector velocity does not have 
a component in the direction of the singular vectors 
associated with the small singular values (notice the 
notch in Fig. 9). This, in fact, is the advantage of hav¬ 
ing the complete SVD available instead of only lim¬ 
ited information on singularities. It should be noted 
that, as a result of passing through a singular configu¬ 
ration. it is possible for the manipulator to switch 
solution branches (for example, going from an elbow- 
up to an elbow-down configuration). 


4.2. A Continuous Version of the Truncated SVD 
Solution 


The use of damped least squares provides the optimal 
solution for tracking a given end-effector trajectory 
under the physical constraints on the joint motions. 
While this solution is optimal, it may be undesirable to 
implement due to the iterative nature of calculating 
the appropriate damping factor A. It can be shown that 
the characteristics of the damped least-squares solu¬ 
tion are very similar to those obtained using the trun¬ 
cated SVD solution. The truncated SVD solution of a 
linear system of equations described by eq. (33), de¬ 
noted here by 0* k \ is defined as 

^=2 jv, (40) 

i-i a i 

where k is an integer less than or equal to the rank r. 
The truncated SVD reduces the solution norm by 
removing all components of the solution that corre¬ 
spond to small singular values while retaining all of 
those associated with larger singular values. The pa¬ 
rameter k is used to define small and large such that a t 
for i k are large, and cr, for i > k are considered 
small. It can be shown that is the minimum resid- 
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Fig. JO. The norm of the 
end-effector velocity tracking 
error , joint velocity , and joint 
accelerations for trajectory B 
using the damped least- 
squares solution at the veloc¬ 
ity level. 


Damped Least Squares 



Trajectory B 


since singular vectors are mutually orthogonal unit 
vectors. 

The advantages of using this form of a solution are 
that, when the SVD is available, it is extremely easy to 
compute. The norms of the truncated singular value 
solutions ||^°||, obtained from eq. (40), are computed 
as i is incremented until either i is equal to the rank or 
the norm is greater than 6^. The latter case is the 
one of interest since it represents reaching the physical 
constraint on the joint velocity. In this case k is now 
known, namely / — 1, and the desired quantity (c — k) 
for which ||0 (c) || is equal to 8^ can be easily obtained 
from (42). This value is then used to obtain the con¬ 
tinuous truncated singular value solution defined by 
(41). An implementation of this algorithm was used to 
simulate control of a PUMA robot for the three trajec¬ 
tories presented above. In all three cases the resultant 
joint trajectories were within 1% of those obtained 
using the damped least-squares solution. A detailed 
error analysis of the difference between the continuous 
truncated SVD solution and the damped least-squares 
solution can be found in Maciejewski (1987). 


ual solution for all 8 in the /c-dimensional subspace 
spanned by Vj for i ^ k (Marquardt 1970). For cases 
where o k » a k+l and X falls between the two largely 
separated singular values o k and cj k + u the results for 
the two types of solutions will be approximately the 
same. By modifying the truncated SVD to be continu¬ 
ous instead of a stepwise function of k , a solution that 
is virtually identical to the damped least-squares solu¬ 
tion can be obtained. This type of solution will be 
denoted 6 (c) and is defined by 

0(0 = y £■ v + (c ~ k)Xk+ ' v (41) 

,rt ffi a k+ , 

where c is a real number less than or equal to the rank, 
and k is the greatest integer less than or equal to c. 

The norm of this type of solution is given by 

l ^ ,||2 -|(f) !+ ( l£ ^r 1 ) 2 ' (42) 


4.3. Resolved Acceleration Control 


The above sections have discussed how to deal with 
hard constraints on the joint angle velocities in the 
presence of singularities by removing those compo¬ 
nents associated with small singular values. In many 
practical cases, however, the joint accelerations will be 
the limiting factor. The same techniques are equally 
applicable, since resolved acceleration control (Luh, 
Walker, and Paul 1980) still requires a solution based 
on some sort of inverse of the Jacobian. This is easily 
seen by differentiating eq. (33) to obtain 

JO +je = \. (43) 

Thus for a given state of the manipulator, the joint 
accelerations required to achieve a desired end-effector 
acceleration can be computed. These joint accelera¬ 
tions can become infinite in the presence of singular¬ 
ities, so that once again a practical solution is to mini¬ 
mize the residual ||x — (J6 + 70)|| under the constraint 
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Fig. 11. The norm of the 
end-effector velocity and ac¬ 
celeration tracking error and 
joint angle velocity and 
acceleration norms for tra¬ 
jectory B using the damped 
least-squares solution for 
resolved acceleration control. 


of physically achievable accelerations defined by ||0|| 
®max- As an example, consider the use of velocity con¬ 
trol for tracking trajectory B given in Figure 10. A plot 
of the joint accelerations required to maintain the 
desired velocities is also given in this figure. The form 
of the joint acceleration norm is typical in that a very 
large spike is located at the singularity, with two 
smaller peaks on either side. Note that this type of 
behavior is not immediately apparent if one considers 
simply differentiating the norm of the joint velocities, 
but it becomes clear after considering the effect of 
passing through a singularity. In particular, the first 
acceleration peak is a result of the joints accelerating 
to match the large velocity required due to the compo¬ 
nent that is in the direction associated with the small 
singular value. The acceleration then goes to zero 
since the hard constraint on the velocities limits the 
effect of this singular component. As the manipulator 
passes through the singularity, however, the compo¬ 
nent along the small singular value switches sign so 
that those joints required to match it must decelerate 
to zero and then accelerate in the opposite direction, 
thus resulting in the large spike in the acceleration 
curve. There is no spike in the velocity curve at this 
point, since the velocity constraint is still in effect, but 
it is now limiting the velocity in the other direction. 
The final peak in the acceleration is then the result of 
leaving the singular region so that the joints deceler¬ 
ate, since large velocities are not required outside of 
singular regions. 

In order to obtain an optimum solution in the pres¬ 
ence of hard constraints on the joint accelerations, the 
continuous truncated SVD technique was applied to 
the solution ofeq. (43). The resulting joint velocities 
and accelerations along with the end-effector tracking 
errors are presented in Figure 11. As expected, the 
end-effector acceleration tracking error is zero, except 
near the singularities, where the hard constraint on the 
acceleration is encountered. This constant limit on the 
acceleration results in the triangular shape of the joint 
velocity curves, which alternate between acceleration 
and deceleration. Imposing such an acceleration con¬ 
straint effectively increases the region in which the 
effect of the singularity is felt. This results in an end- 
effector velocity tracking error that is non-zero for a 
larger portion of the trajectory around the two singular 
configurations. This tracking error, however, is the 
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minimum achievable for the given physical limit on 
the acceleration of the manipulator. 


5. Conclusions 


The conclusions of the work presented here can be 
divided into those relating to calculation of the SVD 
in real time and those pertaining to the advantages of 
having the SVD of the Jacobian available for the con¬ 
trol of robotic manipulators. With regard to the first 
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point, the major advantage of the implementation of 
the SVD algorithm presented here is that by using the 
known SVD of a previous Jacobian, the computa¬ 
tional effort of calculating the SVD of the current Ja¬ 
cobian can be greatly reduced. By applying the rota¬ 
tion implied by the previously known singular vectors, 
a nearly orthogonal matrix results, which will typically 
converge in a single sweep. This fact alone reduces the 
computation time by greater than a factor of three, 
since it has been shown that the original algorithm will 
require between three and four sweeps to converge. In 
addition, the computation time is now known a priori, 
also removing the need for convergence tests. The 
potential difficulty of accumulating the error of pre¬ 
vious computations is removed by alternating the use 
of the input and output singular vectors with which to 
apply the rotations. This modification, while improv¬ 
ing the error characteristics, does not affect the com¬ 
putational requirements. Finally, by taking advantage 
of small angle approximations, the number of floating 
point multiplications required in the plane rotations is 
cut in half, with the total number of floating point 
operations cut by a third. 

The availability of a computationally efficient algo¬ 
rithm for computing the SVD of the Jacobian results 
in a large number of potential applications, including 
the real-time evaluation of dexterity and utilization of 
redundant degrees of freedom. This work illustrates 
the application of the SVD to the fundamental prob¬ 
lem of dealing with singularities. The optimal solution 
to the damped least-squares formulation can be easily 
and efficiently obtained when the SVD is available. 
Furthermore, it has been shown how a continuous 
form of the truncated SVD solution can be used in 
place of the damped least-squares solution. The use of 
this formulation at either the velocity or acceleration 
level allows operation through singular configurations 
without violating physical constraints or generating 
spurious motions. 
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